#function to simulate decomopsition of behavioral response. This one: only use countefactual self-investment policy function
function Simulate_decomp_1(prim_2000::Primitives, prim_2010::Primitives, est::Estimands, res_2000::Results, res_2010::Results, res_2000_cfact::Results, res_2010_cfact::Results, grids::Grids,
    package::Vector{Matrix{Any}}, kid_pctiles::Vector{Float64}; nsims::Int64 = 2000, aopt::Int64 = 0, statcheck::Int64 = 0)
    #@unpack θ, κ_hs, κ_coll, σ_ζ, σ_ε2, ξ, ρ_ha, σ_a, μ_a, ϕ, Δ_1, Δ_2, Δ_3, Δ_4, Δ_5, η_1, η_2, η_3, η_4, σ_η  = est #unpacking
    @unpack θ, κ, σ_ζ, σ_ε2, ξ, ρ_ha, σ_a, μ_a, ϕ, Δ_1, Δ_2, Δ_3, Δ_4, Δ_mod, η_1, η_2, η_3, σ_η, α_1, α_2, α_3, Δ_5  = est #unpacking
    @unpack x_func, t_func, n_func, val_func_2 = res_2000
    @unpack val_func_3 = res_2010_cfact
    @unpack a_grid, ε2_grid, ε3_grid_hs, ε3_grid_coll, h_grid = grids #more unpacking
    @unpack β, skill_prices, prices, nh, na, ne, nS, nl, govt_inv, fips, μ_h, σ_h, s_t_ratios, pops, prox_mat, mig_mat, hc_norm, amenity, divs, regions = prim_2000
    @unpack kid_marr_hs, kid_fert_marr_hs, kid_fert_nonmarr_hs, kid_marr_coll, kid_fert_marr_coll, kid_fert_nonmarr_coll, aid_sched = prim_2000
    @unpack μ_ε3, σ_ε3, taxes, tuitions = prim_2000
    @unpack μ_a_b, μ_a_h, Δ_mod_b, Δ_mod_h, η_1_b, η_1_h, ξ_b, ξ_h = est
    @unpack w_mods, marr_mod_hs, marr_mod_coll, fert_mod_coll_marr, fert_mod_coll_nmarr, fert_mod_hs_marr, fert_mod_hs_nmarr = prim_2000
    @unpack μ_a_r_1, μ_a_r_2, μ_a_r_3, μ_a_r_4 = est


    parent_dist = package[4] #parnet HC distributions and marital probabilities
    ε3_grid = [ε3_grid_hs, ε3_grid_coll]
    marr_collect = [kid_marr_hs, kid_marr_coll]
    fert_marr_collect = [kid_fert_marr_hs, kid_fert_marr_coll]
    fert_nonmarr_collect = [kid_fert_nonmarr_hs, kid_fert_nonmarr_coll]
    kid_marr_dist = package[13]

    amen_1 = amenity
    amen_2 = zeros(50)
    amen_3 = zeros(50)
    if aopt==4 #switch to college gratio
        amen_1 = amen_grid[:, 7]
    elseif aopt>0 && aopt<4
        a_ind = 2*(aopt-1) + 1
        amen_2 = amen_grid[:, a_ind]
        amen_3 = amen_grid[:, a_ind+1]
    end


    #preallocations
    data_simul = Any[] #preallocate
    data_simul_p25 = Any[] #preallocate data for p25 simulations
    unif = Uniform(0,1) #uniform 0 to 1
    dist_ε2 = LogNormal(-σ_ε2^2/2, σ_ε2) #early shock distribution
    dist_ζ = Gumbel(0, σ_ζ) #location preference shock distribution

    #interpolated HC grid and policy functions
    interp_n = interpolate(n_func, (NoInterp(), BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated n function
    interp_n_cfact = interpolate(res_2000_cfact.n_func, (NoInterp(), BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated n function
    interp_x = interpolate(x_func, (BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated x function
    interp_t = interpolate(t_func, (BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated x function
    interp_t_2010 = interpolate(res_2010.t_func, (BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated x function
    interp_t_2010_cfact = interpolate(res_2010_cfact.t_func, (BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated x function
    interp_x_2010 = interpolate(res_2010.x_func, (BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated x function



    if statcheck == 1
        interp_n = interpolate(res_2010.n_func, (NoInterp(), BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated n function
        interp_x = interpolate(res_2010.x_func, (BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated x function
        interp_t = interpolate(res_2010.t_func, (BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated x function
        interp_n_cfact = interpolate(res_2010_cfact.n_func, (NoInterp(), BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated n function
    end


    #parent characterisitcs (education, marriage, etc) by type
    parent_chars = [[1, 1], [1, 2], [1, 3], [1, 4], [2, 1], [2, 2], [2, 3], [2, 4]]

    #loop over states and simulations
    for l = 1:nl
        Random.seed!(l*10)

        row_base = 8*(l-1)

        #probability thresholds to determine parent type
        thresholds = zeros(8)
        thresholds[1] = parent_dist[row_base+1, 3]
        for t = 2:8
            thresholds[t] = thresholds[t-1] + parent_dist[row_base+t, 3]
        end

        region = regions[l]


        #loop over simulations
        for i = 1:nsims

            #draw parent type (college, married)
            draw_p = rand() #between 0 and 1
            type = 1
            for t = 2:7
                if draw_p>thresholds[t-1] && draw_p<thresholds[t]
                    type=t
                end
            end
            if draw_p>thresholds[7]
                type = 8
            end

            con, util = 0, 0 #lifetime consumption and utility

            ####Draw race and associated parametrers
            R_draw = rand()
            R = 1
            if R_draw>parent_dist[row_base+type, 8] && R_draw<(parent_dist[row_base+type, 8] + parent_dist[row_base+type, 9])
                R=2
            elseif R_draw>(parent_dist[row_base+type, 8] + parent_dist[row_base+type, 9])
                R = 3
            end

            #now race is drawn! Unpack other relevant stuff
            #get relevant race stuff
            Δ_mods = [0.0, Δ_mod_b, Δ_mod_h] #check
            μ_a_opts = [μ_a, μ_a_r_2, μ_a_r_3, μ_a_r_4]
            η_opts = [η_1, η_1_b, η_1_h] #check
            Δ_mod_R = Δ_mods[R] #check
            μ_a_R = μ_a_opts[region] #check
            η_R = η_opts[R] #check
            ξ_opts = [ξ, ξ_b, ξ_h]
            ξ_R = ξ_opts[R]

            marr_race_collect = [marr_mod_hs, marr_mod_coll]
            fert_race_marr_collect = [fert_mod_hs_marr, fert_mod_coll_marr]
            fert_race_nmarr_collect = [fert_mod_hs_nmarr, fert_mod_coll_nmarr]


            ###draw parent HC and kid ability###
            h_draw = exp(rand(Normal(parent_dist[row_base+type, 4], parent_dist[row_base+type, 5]))) #draw parent HC
            h = get_index(h_draw, h_grid)
            mean = μ_a_R + ρ_ha*(σ_a/σ_h)*(log(h_draw) - μ_h) #conditional distribution mean
            var = (1-ρ_ha^2)*σ_a^2 #conditional distribution variance
            dist_ap = LogNormal(mean, sqrt(var))
            ap_probs = zeros(na)

            #compute child ability probabilities
            for ap = 1:na
                ap_probs[ap] = pdf(dist_ap, a_grid[ap])
            end
            ap_probs./=sum(ap_probs)
            draw_ap, ap = rand(), 1 #child ability draw

            #determine child abiltiy
            ta1 = ap_probs[1]
            ta2 = ta1 + ap_probs[2]
            ta3 = ta2 + ap_probs[3]
            ta4 = ta3 + ap_probs[4]

            if draw_ap>ta1 && draw_ap<ta2
                ap=2
            elseif draw_ap>ta2 && draw_ap<ta3
                ap=3
            elseif draw_ap>ta3 && draw_ap<ta4
                ap=4
            elseif draw_ap>ta4
                ap=5
            end

            ####compute earnings, investmnet, now with non-assortive HC!###
            S_par = parent_chars[type][1]
            par_marr = parent_chars[type][2]
            m = (par_marr>1)*1
            x_sp, t_sp = 0, 0
            h_sp_draw, h_sp, h_sp_p4 = 0, 0, 0
            earn_sp, earn_sp_p4 = 0, 0

            #adjust if spouse is also working
            if par_marr == 3 || par_marr==4
                h_sp_draw = exp(rand(Normal(parent_dist[row_base+type, 6], parent_dist[row_base+type, 7]))) #draw parent HC
                h_sp_p4 = h_sp_draw * rand(LogNormal(μ_ε3[par_marr-2], σ_ε3[par_marr-2]))
                h_sp = get_index(h_sp_draw, h_grid)
                x_sp = interp_x(h_sp, l, ap, par_marr-2, R)
                t_sp = interp_t(h_sp, l, ap, par_marr-2, R)
                earn_sp =  h_sp_draw* exp(log(skill_prices[l,  par_marr-2]+w_mods[R])) * (1-t_sp)
                earn_sp_p4 = h_sp_p4 * exp(log(skill_prices[l,  par_marr-2]+w_mods[R]))
            end

            if par_marr == 2
                t_sp = 0.13 #spouse time investment if not working -- high!
            end

            x = interp_x(h, l, ap, S_par, R) + x_sp  #parent money investment
            t_ind = interp_t(h, l, ap, S_par, R)  #parent investment
            e3_fam = h_draw*exp(log(skill_prices[l, S_par]+w_mods[R])) * (1-t_ind) + earn_sp  #parent family earnings

            #Need: child ability
            #line_p25 = [type l m S_par h_draw h_sp_draw ap e3_fam]
            #push!(data_simul_p25, line_p25)
            t = t_ind + t_sp #adjust time inputs for kids with married parents

            #compute kid's human capital
            h2_val = ξ_R * a_grid[ap] * (x + (1-ϕ)*govt_inv[l]/prices[l])^(1-ϕ) * (t + ϕ * govt_inv[l]/(prices[l] * s_t_ratios[l] * exp(μ_h + (σ_h^2)/2)))^ϕ #kid human capital
            h2 = get_index(h2_val, h_grid)
            h2_f, h2_c = Int64(floor(h2)), Int64(ceil(h2))
            weight = h2 - h2_f


            #construct kid HS and Copllege exepcted value functions
            kid_val_hs = Base.MathConstants.eulergamma*σ_ζ
            kid_val_coll = Base.MathConstants.eulergamma*σ_ζ

            kid_val_hs_sum = 0
            kid_val_coll_sum = 0

            l2_utils = zeros(nl, nS)

            for lp = 1:nl
                val_lp_hs, val_lp_coll = 0, 0

                #factor in moving cost
                movecost_hs, movecost_coll = 0.0, 0.0
                mig=1
                if lp!=l
                    movecost_hs = 14.327 - Δ_2 * (h2_val) - Δ_3 * prox_mat[l, lp] - Δ_5 * amenity[lp] + Δ_mod_R
                    movecost_coll = 14.327 - Δ_2 * (h2_val) - Δ_3 * prox_mat[l, lp] - Δ_4 - Δ_5 * amenity[lp] + Δ_mod_R
                    mig=2
                end

                val_lp_hs += (1/σ_ζ) * (val_func_2[ap, h2_f, lp, mig, 1, R] * (1-weight) + val_func_2[ap, h2_c, lp, mig, 1, R] * weight)
                val_lp_coll += (1/σ_ζ) * (val_func_2[ap, h2_f, lp, mig, 2, R] * (1-weight) + val_func_2[ap, h2_c, lp, mig, 2, R] * weight)

                val_lp_hs-=(1/σ_ζ)*movecost_hs
                val_lp_coll-=(1/σ_ζ)*movecost_coll

                l2_utils[lp, 1] = val_lp_hs*σ_ζ
                l2_utils[lp, 2] = val_lp_coll*σ_ζ

                kid_val_hs_sum+=exp(val_lp_hs)
                kid_val_coll_sum+=exp(val_lp_coll)
            end
            kid_val_hs += σ_ζ * log(kid_val_hs_sum)
            kid_val_coll += σ_ζ * log(kid_val_coll_sum)

            #COLLEGE DECISION
            #parent HC in period 4
            h4_val = h_draw * rand(LogNormal(μ_ε3[S_par], σ_ε3[S_par]))
            h4 = get_index(h4_val, h_grid)
            fam_earn = h4_val * exp(log(skill_prices[l, S_par]+w_mods[R])) + earn_sp_p4
            grant = Grants(fam_earn, ap, aid_sched)
            draw_college = rand(Normal(0, σ_η)) #college preference draw
            c_hs = h4_val * exp(log(skill_prices[l, S_par]+w_mods[R])) * (1-taxes[l]) / prices[l]
            c_coll = h4_val * exp(log(skill_prices[l, S_par]+w_mods[R])) * (1-taxes[l]) / prices[l] - ((tuitions[l]-grant)/(m*prices[l]))
            S = 1
            threshold_coll = utility(c_hs) + (1+θ)*kid_val_hs - (utility(c_coll) + η_3 * (S_par-1) + (1+θ)*(η_R + η_2 * a_grid[ap] + kid_val_coll))
            if draw_college>threshold_coll
                S = 2
                util += (η_R + η_2*a_grid[ap] + draw_college) #update utility
            end

            ######KID FIRST MOVING DECISION
            loc_shocks = rand(dist_ζ, nl) #draw location shocks
            l2_choice = maximum(l2_utils[:,S].+loc_shocks) #maximizing value
            l2 = findmin(abs.(l2_utils[:,S].+loc_shocks.-l2_choice))[2] #index of choice
            l2 = l #no moving allowd, even though college decision assumed the possibility!
            mig = (l2!=l) + 1
            n = interp_n_cfact(ap, h2, l2, mig, S, R) #kid's self-investment choice
            time_coll = (2/9) * (S-1)
            e2p = h2_val*exp(log(skill_prices[l2, S]+w_mods[R]))*(1-n - time_coll) #early life earnings

            ###Determine Moving Decision###
            con += (e2p * (1-taxes[l2]))/prices[l2]
            util += loc_shocks[l2] #utilty shock from location choice
            util -= (Δ_1 - Δ_2 * (h2_val) - Δ_3 * prox_mat[l, l2] - Δ_4*(S-1) - Δ_5 * amenity[l2]) * (l2!=l) #moving cost if move
            util += utility((e2p * (1-taxes[l2]))/prices[l2]) + (α_1 * amen_1[l2] + α_2 * amen_2[l2] + α_3 * amen_3[l2])

            ###Determine Moving Decision###

            #actually draw next-period HC
            h3p_val = rand(dist_ε2)*(a_grid[ap]*(n*h2_val)^κ + h2_val)
            h3p = get_index(h3p_val, h_grid)
            h3_ind_f, h3_ind_c = Int64(floor(h3p)), Int64(ceil(h3p))
            weight = h3p - h3_ind_f

            kid_marr = marr_collect[S]
            kid_fert_marr = fert_marr_collect[S]
            kid_fert_nonmarr = fert_nonmarr_collect[S]
            marr_R_mod = marr_race_collect[S][R]
            fert_marr_R_mod = fert_race_marr_collect[S][R]
            fert_nmarr_R_mod = fert_race_nmarr_collect[S][R]


            #compute marriage/fertility probabilities
            prob_marr = 𝚽(kid_marr[l2, 2] + h3p_val*kid_marr[l2, 3] +  h3p_val^2*kid_marr[l2, 4] +  h3p_val^3*kid_marr[l2, 5] + kid_marr[l2, 7] * (l!=l2) + marr_R_mod)
            prob_fert_marr = 𝚽(kid_fert_marr[l2,2] + h3p_val*kid_fert_marr[l2,3] +  h3p_val^2*kid_fert_marr[l2,4] +  h3p_val^3*kid_fert_marr[l2,5] + fert_marr_R_mod)
            prob_fert_nonmarr = 𝚽(kid_fert_nonmarr[l2,2] + h3p_val*kid_fert_nonmarr[l2,3] +  h3p_val^2*kid_fert_nonmarr[l2,4] +  h3p_val^3*kid_fert_nonmarr[l2,5] + fert_nmarr_R_mod)

            #cap
            if h3p_val>4
                prob_marr = kid_marr[l2, 6]
                prob_fert_marr, prob_fert_nonmarr = kid_fert_marr[l2, 6], kid_fert_nonmarr[l2, 6]
            end

            #get probabilities of grandkids ability
            app_probs = zeros(na, 4)
            for reg = 1:4
                μ_a_R_p = μ_a_opts[reg] #check
                mean = μ_a_R_p + ρ_ha*(σ_a/σ_h)*(log(h3p_val) - μ_h) #conditional distribution mean
                var = (1-ρ_ha^2)*σ_a^2 #conditional distribution variance
                dist_app = LogNormal(mean, sqrt(var))

                #compute child ability probabilities
                for app = 1:na
                    app_probs[app, reg] = pdf(dist_app, a_grid[app])
                end
                app_probs[:, reg]./=sum(app_probs[:, reg])
            end

            l_utils = zeros(nl) #compute expected utility from each possible location
            l_utils_norec = zeros(nl)


            for lp = 1:nl

                regp = regions[lp]
                val_lp = 0.0
                val_lp_norec = 0.0

                #factor in moving cost
                movecost = 0.0
                if lp!=l2
                    #movecost = Δ_1 - Δ_2 * (h3p_val-hc_norm) - Δ_3 * prox_mat[l, lp] - Δ_4 * (S-1) + Δ_5 * pops[lp]
                    #movecost = Δ_1 - Δ_2 * (h3p_val) - Δ_3 * prox_mat[l, lp] - Δ_4 * (S-1) + Δ_mod #adjusted for higher period-3 moving cost
                    movecost = 750 - Δ_2 * (h3p_val) - Δ_3 * prox_mat[l2, lp] - Δ_4 * (S-1) - Δ_5 * amenity[lp] + Δ_mod_R #adjusted for higher period-3 moving cost
                end

                for app = 1:na

                    #####add stuff to continuation value
                    #kids, non-married
                    val_lp += val_func_3[h3_ind_f, lp, app, 1, S, R] * app_probs[app, regp] * (1-prob_marr) * prob_fert_nonmarr * (1-weight)
                    val_lp += val_func_3[h3_ind_c, lp, app, 1, S, R] * app_probs[app, regp] * (1-prob_marr) * prob_fert_nonmarr * (weight)

                    #no kid, no marriage
                    val_lp += val_func_3[h3_ind_f, lp, app, 2, S, R] * app_probs[app, regp] * (1-prob_marr) * (1-prob_fert_nonmarr) * (1-weight)
                    val_lp += val_func_3[h3_ind_c, lp, app, 2, S, R] * app_probs[app, regp] * (1-prob_marr) * (1-prob_fert_nonmarr) * (weight)

                    #kids, married
                    val_lp += val_func_3[h3_ind_f, lp, app, 1, S, R] * app_probs[app, regp] * (prob_marr) * prob_fert_marr * (1-weight) * (1+θ)
                    val_lp += val_func_3[h3_ind_c, lp, app, 1, S, R] * app_probs[app, regp] * (prob_marr) * prob_fert_marr * (weight) * (1+θ)

                    #no kid, marriage
                    val_lp += val_func_3[h3_ind_f, lp, app, 2, S, R] * app_probs[app, regp] * (prob_marr) * (1-prob_fert_marr) * (1-weight) * (1+θ)
                    val_lp += val_func_3[h3_ind_c, lp, app, 2, S, R] * app_probs[app, regp] * (prob_marr) * (1-prob_fert_marr) * (weight) * (1+θ)


                    ###############

                    #val_lp_norec += res_2000.val_func_3[h3_ind_f, lp, app, 1, S] * app_probs[app] * (1-prob_marr) * prob_fert_nonmarr * (1-weight)
                    #val_lp_norec += res_2000.val_func_3[h3_ind_c, lp, app, 1, S] * app_probs[app] * (1-prob_marr) * prob_fert_nonmarr * (weight)

                    #no kid, no marriage
                    #val_lp_norec += res_2000.val_func_3[h3_ind_f, lp, app, 2, S] * app_probs[app] * (1-prob_marr) * (1-prob_fert_nonmarr) * (1-weight)
                    #val_lp_norec += res_2000.val_func_3[h3_ind_c, lp, app, 2, S] * app_probs[app] * (1-prob_marr) * (1-prob_fert_nonmarr) * (weight)

                    #kids, married
                    #val_lp_norec += res_2000.val_func_3[h3_ind_f, lp, app, 1, S] * app_probs[app] * (prob_marr) * prob_fert_marr * (1-weight) * (1+θ)
                    #val_lp_norec += res_2000.val_func_3[h3_ind_c, lp, app, 1, S] * app_probs[app] * (prob_marr) * prob_fert_marr * (weight) * (1+θ)

                    #no kid, marriage
                    #val_lp_norec += res_2000.val_func_3[h3_ind_f, lp, app, 2, S] * app_probs[app] * (prob_marr) * (1-prob_fert_marr) * (1-weight) * (1+θ)
                    #val_lp_norec += res_2000.val_func_3[h3_ind_c, lp, app, 2, S] * app_probs[app] * (prob_marr) * (1-prob_fert_marr) * (weight) * (1+θ)

                end
                val_lp -= movecost
                val_lp_norec -= movecost
                l_utils[lp] = val_lp
                l_utils_norec[lp] = val_lp_norec
            end

            loc_shocks = rand(dist_ζ, nl) #draw location shocks
            lp_choice = maximum(l_utils.+loc_shocks) #maximizing value
            lp = findmin(abs.(l_utils.+loc_shocks.-lp_choice))[2] #index of choice
            region_p = regions[lp]
            μ_a_R_p = μ_a_opts[region_p] #check

            lp_choice_norec = maximum(l_utils_norec.+loc_shocks) #maximizing value
            lp_norec = findmin(abs.(l_utils_norec.+loc_shocks.-lp_choice_norec))[2] #index of choice

            #update utility
            util += β * loc_shocks[lp] #location shock
            util -= (Δ_mod - Δ_2 * (h3p_val) - Δ_3 * prox_mat[l2, lp] - Δ_4 * (S-1) - Δ_5 * amenity[lp]) * (l2!=lp) * β #moving cost

            #####draw marriage/fertility outcomes
            ##marriag
            #draw
            marr_draw = rand()
            mp = 1 #non-married
            if marr_draw<prob_marr
                mp=2 #married
            end

            ##fertility
            prob_fert = prob_fert_marr
            if mp == 1 #non-married
                prob_fert = prob_fert_nonmarr
            end

            fert_draw = rand()
            fert = 0 #no kids
            if fert_draw<prob_fert
                fert = 1 #kids
            end

            #draw grandchild ability
            draw_app, app = rand(), 1 #child ability draw

            #determine child abiltiy
            tap1 = app_probs[1, region_p]
            tap2 = tap1 + app_probs[2, region_p]
            tap3 = tap2 + app_probs[3, region_p]
            tap4 = tap3 + app_probs[4, region_p]

            if draw_app>tap1 && draw_app<tap2
                app=2
            elseif draw_app>tap2 && draw_app<tap3
                app=3
            elseif draw_app>tap3 && draw_app<tap4
                app=4
            elseif draw_app>tap4
                app=5
            end

            tp = interp_t_2010_cfact(h3p, lp, app, S, R) #parent investment
            xp = interp_x_2010(h3p, lp, app, S, R) #parent goods


            #compute earnings with more flexible child marriage setup. Could technically add stuff to this, but this is already pretty grueling
            kid_marr_draw = rand()
            row_kid_marr_dist = 6*(l2-1) + 3*(S-1)
            spouse_kid_type, μ_hc_sp, σ_hc_sp = 1, 0, 0
            threshold_kid_marr_1 = kid_marr_dist[row_kid_marr_dist+1, 4]
            threshold_kid_marr_2 = threshold_kid_marr_1 + kid_marr_dist[row_kid_marr_dist+2, 4]
            if kid_marr_draw>threshold_kid_marr_1 && kid_marr_draw<threshold_kid_marr_2
                spouse_kid_type, μ_hc_sp, σ_hc_sp = 2, kid_marr_dist[row_kid_marr_dist+2, 5], kid_marr_dist[row_kid_marr_dist+2, 6]
            elseif kid_marr_draw>threshold_kid_marr_2
                spouse_kid_type, μ_hc_sp, σ_hc_sp = 3, kid_marr_dist[row_kid_marr_dist+3, 5], kid_marr_dist[row_kid_marr_dist+3, 6]
            end

            earn_sp_kid = 0
            earn_sp_kid_p4 = 0
            x_sp_kid = 0
            time_sp_kid = 0.13
            hc_draw_sp = rand(LogNormal(μ_hc_sp, σ_hc_sp))
            h3p_sp = get_index(hc_draw_sp, h_grid)
            h3p_shocks_sp = [rand(LogNormal(μ_ε3[1], σ_ε3[1])), rand(LogNormal(μ_ε3[2], σ_ε3[2]))]
            if spouse_kid_type>1 #update spousal edarnings if working
                time_sp_kid = interp_t_2010_cfact(h3p_sp, lp, app, spouse_kid_type-1, R)
                x_sp_kid = interp_x_2010(h3p_sp, lp, app, spouse_kid_type-1, R)
                earn_sp_kid = hc_draw_sp * exp(log(prim_2010.skill_prices[lp, spouse_kid_type-1]+w_mods[R])) * (1-time_sp_kid * fert)
                earn_sp_kid_p4 = exp(log(skill_prices[lp, spouse_kid_type-1]+w_mods[R])) * h3p_sp * h3p_shocks_sp[spouse_kid_type-1]
            end

            skill_price_starting = exp(log(prim_2010.skill_prices[l, S]+w_mods[R]))
            skill_price_ending = exp(log(prim_2010.skill_prices[lp, S]+w_mods[R]))
            e3p_fam = h3p_val*skill_price_ending * (1-tp*fert) + earn_sp_kid * (mp-1) #child family earnings
            e3p = h3p_val*skill_price_ending * (1-tp*fert) #compute next-period individual earnings

            ########udpate utility and consumption again -- this time a bit involved . . .#####
            #period 3 consumption, utility
            ctemp = (e3p * (1-taxes[lp]))/prices[lp] - xp
            if fert == 1
                ctemp/=1.5 #equivalence scale adjustment
            end
            con += ctemp
            util += (utility(ctemp) + (α_1 * amen_1[lp] + α_2 * amen_2[lp] + α_3 * amen_3[lp])) * (1 + θ*(mp-1)) * β #add to utility, with altruistic adjustment if married

            #realize period-4 HC
            h4p_val = h3p_val * rand(LogNormal(μ_ε3[S], σ_ε3[S]))
            h4p = get_index(h4p_val, h_grid)
            draw_college_g = rand(Normal(0, σ_η)) #college preference draw for grandkid


            e3p_fam = h3p_val*prim_2010.skill_prices[lp, S] * (1-tp*fert) + earn_sp_kid * (mp-1) #child family earnings
            #e3p_fam = h3p_val*prim_2010.skill_prices[lp, S] * (1-tp*fert) * mp #child family earnings
            e3p = h3p_val*prim_2010.skill_prices[lp, S] * (1-tp*fert) #compute next-period individual earnings
            #e3p = e3p_fam/mp #compute next-period individual earnings
            mig = (l!=lp) * 1.0

            #additional migration variables
            mig2 = (l2!=lp) * 1.0 #migrated in period 2-3

            ###update simulated data vector
            line = [fips[l] fips[lp] prim_2010.skill_prices[l, S] prim_2010.skill_prices[lp, S] m mp t e2p e3p h3p_val ap S_par S pops[l] mig type mig2 fips[lp_norec] e3_fam e3p_fam]
            push!(data_simul, line)
        end
    end

    #wrap up and export
    data_simul = vcat(data_simul...) #stick all together into one array
    data_simul_p25 = vcat(data_simul_p25...) #stick all together into one array
    pop_weights = fweights(data_simul[:,14]) #define population weights

    #compute kid and parent percentiles, with these being the final two entries of the data
    ncol = length(data_simul[1,:])

    nrows = length(data_simul[:,1]) #useful for later
    ranks = zeros(nrows, 1)
    for i = 1:nrows
        ranks[i,1] = findmin(abs.(kid_pctiles .- data_simul[i,ncol]))[2]
    end

    data_simul = hcat(data_simul, ranks)
    #data_simul_p25 = hcat(data_simul_p25, ranks[:,1])
    output_simul = data_simul
    output_simul #return simulated output
end




#function to simulate decomopsition of behavioral response. This one: only use countefactual self-investment policy function
function Simulate_decomp_2(prim_2000::Primitives, prim_2010::Primitives, est::Estimands, res_2000::Results, res_2010::Results, res_2000_cfact::Results, res_2010_cfact::Results, grids::Grids,
    package::Vector{Matrix{Any}}, kid_pctiles::Vector{Float64}; nsims::Int64 = 2000, aopt::Int64 = 0, statcheck::Int64 = 0)
    #@unpack θ, κ_hs, κ_coll, σ_ζ, σ_ε2, ξ, ρ_ha, σ_a, μ_a, ϕ, Δ_1, Δ_2, Δ_3, Δ_4, Δ_5, η_1, η_2, η_3, η_4, σ_η  = est #unpacking
    @unpack θ, κ, σ_ζ, σ_ε2, ξ, ρ_ha, σ_a, μ_a, ϕ, Δ_1, Δ_2, Δ_3, Δ_4, Δ_mod, η_1, η_2, η_3, σ_η, α_1, α_2, α_3, Δ_5  = est #unpacking
    @unpack x_func, t_func, n_func = res_2000
    @unpack val_func_3 = res_2010_cfact
    @unpack val_func_2 = res_2000_cfact
    @unpack a_grid, ε2_grid, ε3_grid_hs, ε3_grid_coll, h_grid = grids #more unpacking
    @unpack β, skill_prices, prices, nh, na, ne, nS, nl, govt_inv, fips, μ_h, σ_h, s_t_ratios, pops, prox_mat, mig_mat, hc_norm, amenity, divs, regions = prim_2000
    @unpack kid_marr_hs, kid_fert_marr_hs, kid_fert_nonmarr_hs, kid_marr_coll, kid_fert_marr_coll, kid_fert_nonmarr_coll, aid_sched = prim_2000
    @unpack μ_ε3, σ_ε3, taxes, tuitions = prim_2000
    @unpack μ_a_b, μ_a_h, Δ_mod_b, Δ_mod_h, η_1_b, η_1_h, ξ_b, ξ_h = est
    @unpack w_mods, marr_mod_hs, marr_mod_coll, fert_mod_coll_marr, fert_mod_coll_nmarr, fert_mod_hs_marr, fert_mod_hs_nmarr = prim_2000
    @unpack μ_a_r_1, μ_a_r_2, μ_a_r_3, μ_a_r_4 = est


    parent_dist = package[4] #parnet HC distributions and marital probabilities
    ε3_grid = [ε3_grid_hs, ε3_grid_coll]
    marr_collect = [kid_marr_hs, kid_marr_coll]
    fert_marr_collect = [kid_fert_marr_hs, kid_fert_marr_coll]
    fert_nonmarr_collect = [kid_fert_nonmarr_hs, kid_fert_nonmarr_coll]
    kid_marr_dist = package[13]

    #preallocations
    data_simul = Any[] #preallocate
    data_simul_p25 = Any[] #preallocate data for p25 simulations
    unif = Uniform(0,1) #uniform 0 to 1
    dist_ε2 = LogNormal(-σ_ε2^2/2, σ_ε2) #early shock distribution
    dist_ζ = Gumbel(0, σ_ζ) #location preference shock distribution

    amen_1 = amenity
    amen_2 = zeros(50)
    amen_3 = zeros(50)
    if aopt==4 #switch to college gratio
        amen_1 = amen_grid[:, 7]
    elseif aopt>0 && aopt<4
        a_ind = 2*(aopt-1) + 1
        amen_2 = amen_grid[:, a_ind]
        amen_3 = amen_grid[:, a_ind+1]
    end

    #interpolated HC grid and policy functions
    interp_n = interpolate(n_func, (NoInterp(), BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated n function
    interp_n_cfact = interpolate(res_2000_cfact.n_func, (NoInterp(), BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated n function
    interp_x = interpolate(x_func, (BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated x function
    interp_t = interpolate(t_func, (BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated x function
    interp_t_2010 = interpolate(res_2010.t_func, (BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated x function
    interp_t_2010_cfact = interpolate(res_2010_cfact.t_func, (BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated x function
    interp_x_2010 = interpolate(res_2010.x_func, (BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated x function

    if statcheck == 1
        interp_n = interpolate(res_2010.n_func, (NoInterp(), BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated n function
        interp_x = interpolate(res_2010.x_func, (BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated x function
        interp_t = interpolate(res_2010.t_func, (BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated x function
        interp_n_cfact = interpolate(res_2010_cfact.n_func, (NoInterp(), BSpline(Linear()), NoInterp(), NoInterp(), NoInterp(), NoInterp())) #interpolated n function
    end

    #parent characterisitcs (education, marriage, etc) by type
    parent_chars = [[1, 1], [1, 2], [1, 3], [1, 4], [2, 1], [2, 2], [2, 3], [2, 4]]

    #loop over states and simulations
    for l = 1:nl
        Random.seed!(l*10)

        row_base = 8*(l-1)

        #probability thresholds to determine parent type
        thresholds = zeros(8)
        thresholds[1] = parent_dist[row_base+1, 3]
        for t = 2:8
            thresholds[t] = thresholds[t-1] + parent_dist[row_base+t, 3]
        end

        region = regions[l]


        #loop over simulations
        for i = 1:nsims


            #draw parent type (college, married)
            draw_p = rand() #between 0 and 1
            type = 1
            for t = 2:7
                if draw_p>thresholds[t-1] && draw_p<thresholds[t]
                    type=t
                end
            end
            if draw_p>thresholds[7]
                type = 8
            end

            con, util = 0, 0 #lifetime consumption and utility

            ####Draw race and associated parametrers
            R_draw = rand()
            R = 1
            if R_draw>parent_dist[row_base+type, 8] && R_draw<(parent_dist[row_base+type, 8] + parent_dist[row_base+type, 9])
                R=2
            elseif R_draw>(parent_dist[row_base+type, 8] + parent_dist[row_base+type, 9])
                R = 3
            end

            #now race is drawn! Unpack other relevant stuff
            #get relevant race stuff
            Δ_mods = [0.0, Δ_mod_b, Δ_mod_h] #check
            μ_a_opts = [μ_a, μ_a_r_2, μ_a_r_3, μ_a_r_4]
            η_opts = [η_1, η_1_b, η_1_h] #check
            Δ_mod_R = Δ_mods[R] #check
            μ_a_R = μ_a_opts[region] #check
            η_R = η_opts[R] #check
            ξ_opts = [ξ, ξ_b, ξ_h]
            ξ_R = ξ_opts[R]

            marr_race_collect = [marr_mod_hs, marr_mod_coll]
            fert_race_marr_collect = [fert_mod_hs_marr, fert_mod_coll_marr]
            fert_race_nmarr_collect = [fert_mod_hs_nmarr, fert_mod_coll_nmarr]


            ###draw parent HC and kid ability###
            h_draw = exp(rand(Normal(parent_dist[row_base+type, 4], parent_dist[row_base+type, 5]))) #draw parent HC
            h = get_index(h_draw, h_grid)
            mean = μ_a_R + ρ_ha*(σ_a/σ_h)*(log(h_draw) - μ_h) #conditional distribution mean
            var = (1-ρ_ha^2)*σ_a^2 #conditional distribution variance
            dist_ap = LogNormal(mean, sqrt(var))
            ap_probs = zeros(na)

            #compute child ability probabilities
            for ap = 1:na
                ap_probs[ap] = pdf(dist_ap, a_grid[ap])
            end
            ap_probs./=sum(ap_probs)
            draw_ap, ap = rand(), 1 #child ability draw

            #determine child abiltiy
            ta1 = ap_probs[1]
            ta2 = ta1 + ap_probs[2]
            ta3 = ta2 + ap_probs[3]
            ta4 = ta3 + ap_probs[4]

            if draw_ap>ta1 && draw_ap<ta2
                ap=2
            elseif draw_ap>ta2 && draw_ap<ta3
                ap=3
            elseif draw_ap>ta3 && draw_ap<ta4
                ap=4
            elseif draw_ap>ta4
                ap=5
            end

            ####compute earnings, investmnet, now with non-assortive HC!###
            S_par = parent_chars[type][1]
            par_marr = parent_chars[type][2]
            m = (par_marr>1)*1
            x_sp, t_sp = 0, 0
            h_sp_draw, h_sp, h_sp_p4 = 0, 0, 0
            earn_sp, earn_sp_p4 = 0, 0

            #adjust if spouse is also working
            if par_marr == 3 || par_marr==4
                h_sp_draw = exp(rand(Normal(parent_dist[row_base+type, 6], parent_dist[row_base+type, 7]))) #draw parent HC
                h_sp_p4 = h_sp_draw * rand(LogNormal(μ_ε3[par_marr-2], σ_ε3[par_marr-2]))
                h_sp = get_index(h_sp_draw, h_grid)
                x_sp = interp_x(h_sp, l, ap, par_marr-2, R)
                t_sp = interp_t(h_sp, l, ap, par_marr-2, R)
                earn_sp =  h_sp_draw* exp(log(skill_prices[l,  par_marr-2]+w_mods[R])) * (1-t_sp)
                earn_sp_p4 = h_sp_p4 * exp(log(skill_prices[l,  par_marr-2]+w_mods[R]))
            end

            if par_marr == 2
                t_sp = 0.13 #spouse time investment if not working -- high!
            end

            x = interp_x(h, l, ap, S_par, R) + x_sp  #parent money investment
            t_ind = interp_t(h, l, ap, S_par, R)  #parent investment
            e3_fam = h_draw*exp(log(skill_prices[l, S_par]+w_mods[R])) * (1-t_ind) + earn_sp  #parent family earnings

            #Need: child ability
            #line_p25 = [type l m S_par h_draw h_sp_draw ap e3_fam]
            #push!(data_simul_p25, line_p25)
            t = t_ind + t_sp #adjust time inputs for kids with married parents

            #compute kid's human capital
            h2_val = ξ_R * a_grid[ap] * (x + (1-ϕ)*govt_inv[l]/prices[l])^(1-ϕ) * (t + ϕ * govt_inv[l]/(prices[l] * s_t_ratios[l] * exp(μ_h + (σ_h^2)/2)))^ϕ #kid human capital
            h2 = get_index(h2_val, h_grid)
            h2_f, h2_c = Int64(floor(h2)), Int64(ceil(h2))
            weight = h2 - h2_f


            #construct kid HS and Copllege exepcted value functions
            kid_val_hs = Base.MathConstants.eulergamma*σ_ζ
            kid_val_coll = Base.MathConstants.eulergamma*σ_ζ

            kid_val_hs_sum = 0
            kid_val_coll_sum = 0

            l2_utils = zeros(nl, nS)

            for lp = 1:nl
                val_lp_hs, val_lp_coll = 0, 0

                #factor in moving cost
                movecost_hs, movecost_coll = 0.0, 0.0
                mig=1
                if lp!=l
                    movecost_hs = 750 - Δ_2 * (h2_val) - Δ_3 * prox_mat[l, lp] - Δ_5 * amenity[lp] + Δ_mod_R
                    movecost_coll = 750 - Δ_2 * (h2_val) - Δ_3 * prox_mat[l, lp] - Δ_4 - Δ_5 * amenity[lp] + Δ_mod_R
                    mig=2
                end

                val_lp_hs += (1/σ_ζ) * (val_func_2[ap, h2_f, lp, mig, 1, R] * (1-weight) + val_func_2[ap, h2_c, lp, mig, 1, R] * weight)
                val_lp_coll += (1/σ_ζ) * (val_func_2[ap, h2_f, lp, mig, 2, R] * (1-weight) + val_func_2[ap, h2_c, lp, mig, 2, R] * weight)

                val_lp_hs-=(1/σ_ζ)*movecost_hs
                val_lp_coll-=(1/σ_ζ)*movecost_coll

                l2_utils[lp, 1] = val_lp_hs*σ_ζ
                l2_utils[lp, 2] = val_lp_coll*σ_ζ

                kid_val_hs_sum+=exp(val_lp_hs)
                kid_val_coll_sum+=exp(val_lp_coll)
            end
            kid_val_hs += σ_ζ * log(kid_val_hs_sum)
            kid_val_coll += σ_ζ * log(kid_val_coll_sum)

            #COLLEGE DECISION
            #parent HC in period 4
            h4_val = h_draw * rand(LogNormal(μ_ε3[S_par], σ_ε3[S_par]))
            h4 = get_index(h4_val, h_grid)
            fam_earn = h4_val * exp(log(skill_prices[l, S_par]+w_mods[R])) + earn_sp_p4
            grant = Grants(fam_earn, ap, aid_sched)
            draw_college = rand(Normal(0, σ_η)) #college preference draw
            c_hs = h4_val * exp(log(skill_prices[l, S_par]+w_mods[R])) * (1-taxes[l]) / prices[l]
            c_coll = h4_val * exp(log(skill_prices[l, S_par]+w_mods[R])) * (1-taxes[l]) / prices[l] - ((tuitions[l]-grant)/(m*prices[l]))
            S = 1
            threshold_coll = utility(c_hs) + (1+θ)*kid_val_hs - (utility(c_coll) + η_3 * (S_par-1) + (1+θ)*(η_R + η_2 * a_grid[ap] + kid_val_coll))
            if draw_college>threshold_coll
                S = 2
                util += (η_R + η_2*a_grid[ap] + draw_college) #update utility
            end

            ######KID FIRST MOVING DECISION
            loc_shocks = rand(dist_ζ, nl) #draw location shocks
            l2_choice = maximum(l2_utils[:,S].+loc_shocks) #maximizing value
            l2 = findmin(abs.(l2_utils[:,S].+loc_shocks.-l2_choice))[2] #index of choice
            l2 = l #no moving allowd, even though college decision assumed the possibility!
            mig = (l2!=l) + 1
            n = interp_n_cfact(ap, h2, l2, mig, S, R) #kid's self-investment choice
            time_coll = (2/9) * (S-1)
            e2p = h2_val*exp(log(skill_prices[l2, S]+w_mods[R]))*(1-n - time_coll) #early life earnings

            ###Determine Moving Decision###
            con += (e2p * (1-taxes[l2]))/prices[l2]
            util += loc_shocks[l2] #utilty shock from location choice
            util -= (Δ_1 - Δ_2 * (h2_val) - Δ_3 * prox_mat[l, l2] - Δ_4*(S-1) - Δ_5 * amenity[l2]) * (l2!=l) #moving cost if move
            util += utility((e2p * (1-taxes[l2]))/prices[l2]) + (α_1 * amen_1[l2] + α_2 * amen_2[l2] + α_3 * amen_3[l2])

            ###Determine Moving Decision###

            #actually draw next-period HC
            h3p_val = rand(dist_ε2)*(a_grid[ap]*(n*h2_val)^κ + h2_val)
            h3p = get_index(h3p_val, h_grid)
            h3_ind_f, h3_ind_c = Int64(floor(h3p)), Int64(ceil(h3p))
            weight = h3p - h3_ind_f

            kid_marr = marr_collect[S]
            kid_fert_marr = fert_marr_collect[S]
            kid_fert_nonmarr = fert_nonmarr_collect[S]
            marr_R_mod = marr_race_collect[S][R]
            fert_marr_R_mod = fert_race_marr_collect[S][R]
            fert_nmarr_R_mod = fert_race_nmarr_collect[S][R]


            #compute marriage/fertility probabilities
            prob_marr = 𝚽(kid_marr[l2, 2] + h3p_val*kid_marr[l2, 3] +  h3p_val^2*kid_marr[l2, 4] +  h3p_val^3*kid_marr[l2, 5] + kid_marr[l2, 7] * (l!=l2) + marr_R_mod)
            prob_fert_marr = 𝚽(kid_fert_marr[l2,2] + h3p_val*kid_fert_marr[l2,3] +  h3p_val^2*kid_fert_marr[l2,4] +  h3p_val^3*kid_fert_marr[l2,5] + fert_marr_R_mod)
            prob_fert_nonmarr = 𝚽(kid_fert_nonmarr[l2,2] + h3p_val*kid_fert_nonmarr[l2,3] +  h3p_val^2*kid_fert_nonmarr[l2,4] +  h3p_val^3*kid_fert_nonmarr[l2,5] + fert_nmarr_R_mod)

            #cap
            if h3p_val>4
                prob_marr = kid_marr[l2, 6]
                prob_fert_marr, prob_fert_nonmarr = kid_fert_marr[l2, 6], kid_fert_nonmarr[l2, 6]
            end

            #get probabilities of grandkids ability
            app_probs = zeros(na, 4)
            for reg = 1:4
                μ_a_R_p = μ_a_opts[reg] #check
                mean = μ_a_R_p + ρ_ha*(σ_a/σ_h)*(log(h3p_val) - μ_h) #conditional distribution mean
                var = (1-ρ_ha^2)*σ_a^2 #conditional distribution variance
                dist_app = LogNormal(mean, sqrt(var))

                #compute child ability probabilities
                for app = 1:na
                    app_probs[app, reg] = pdf(dist_app, a_grid[app])
                end
                app_probs[:, reg]./=sum(app_probs[:, reg])
            end

            l_utils = zeros(nl) #compute expected utility from each possible location
            l_utils_norec = zeros(nl)


            for lp = 1:nl

                regp = regions[lp]
                val_lp = 0.0
                val_lp_norec = 0.0

                #factor in moving cost
                movecost = 0.0
                if lp!=l2
                    #movecost = Δ_1 - Δ_2 * (h3p_val-hc_norm) - Δ_3 * prox_mat[l, lp] - Δ_4 * (S-1) + Δ_5 * pops[lp]
                    #movecost = Δ_1 - Δ_2 * (h3p_val) - Δ_3 * prox_mat[l, lp] - Δ_4 * (S-1) + Δ_mod #adjusted for higher period-3 moving cost
                    movecost = 750 - Δ_2 * (h3p_val) - Δ_3 * prox_mat[l2, lp] - Δ_4 * (S-1) - Δ_5 * amenity[lp] + Δ_mod_R #adjusted for higher period-3 moving cost
                end

                for app = 1:na

                    #####add stuff to continuation value
                    #kids, non-married
                    val_lp += val_func_3[h3_ind_f, lp, app, 1, S, R] * app_probs[app, regp] * (1-prob_marr) * prob_fert_nonmarr * (1-weight)
                    val_lp += val_func_3[h3_ind_c, lp, app, 1, S, R] * app_probs[app, regp] * (1-prob_marr) * prob_fert_nonmarr * (weight)

                    #no kid, no marriage
                    val_lp += val_func_3[h3_ind_f, lp, app, 2, S, R] * app_probs[app, regp] * (1-prob_marr) * (1-prob_fert_nonmarr) * (1-weight)
                    val_lp += val_func_3[h3_ind_c, lp, app, 2, S, R] * app_probs[app, regp] * (1-prob_marr) * (1-prob_fert_nonmarr) * (weight)

                    #kids, married
                    val_lp += val_func_3[h3_ind_f, lp, app, 1, S, R] * app_probs[app, regp] * (prob_marr) * prob_fert_marr * (1-weight) * (1+θ)
                    val_lp += val_func_3[h3_ind_c, lp, app, 1, S, R] * app_probs[app, regp] * (prob_marr) * prob_fert_marr * (weight) * (1+θ)

                    #no kid, marriage
                    val_lp += val_func_3[h3_ind_f, lp, app, 2, S, R] * app_probs[app, regp] * (prob_marr) * (1-prob_fert_marr) * (1-weight) * (1+θ)
                    val_lp += val_func_3[h3_ind_c, lp, app, 2, S, R] * app_probs[app, regp] * (prob_marr) * (1-prob_fert_marr) * (weight) * (1+θ)


                    ###############

                    #val_lp_norec += res_2000.val_func_3[h3_ind_f, lp, app, 1, S] * app_probs[app] * (1-prob_marr) * prob_fert_nonmarr * (1-weight)
                    #val_lp_norec += res_2000.val_func_3[h3_ind_c, lp, app, 1, S] * app_probs[app] * (1-prob_marr) * prob_fert_nonmarr * (weight)

                    #no kid, no marriage
                    #val_lp_norec += res_2000.val_func_3[h3_ind_f, lp, app, 2, S] * app_probs[app] * (1-prob_marr) * (1-prob_fert_nonmarr) * (1-weight)
                    #val_lp_norec += res_2000.val_func_3[h3_ind_c, lp, app, 2, S] * app_probs[app] * (1-prob_marr) * (1-prob_fert_nonmarr) * (weight)

                    #kids, married
                    #val_lp_norec += res_2000.val_func_3[h3_ind_f, lp, app, 1, S] * app_probs[app] * (prob_marr) * prob_fert_marr * (1-weight) * (1+θ)
                    #val_lp_norec += res_2000.val_func_3[h3_ind_c, lp, app, 1, S] * app_probs[app] * (prob_marr) * prob_fert_marr * (weight) * (1+θ)

                    #no kid, marriage
                    #val_lp_norec += res_2000.val_func_3[h3_ind_f, lp, app, 2, S] * app_probs[app] * (prob_marr) * (1-prob_fert_marr) * (1-weight) * (1+θ)
                    #val_lp_norec += res_2000.val_func_3[h3_ind_c, lp, app, 2, S] * app_probs[app] * (prob_marr) * (1-prob_fert_marr) * (weight) * (1+θ)

                end
                val_lp -= movecost
                val_lp_norec -= movecost
                l_utils[lp] = val_lp
                l_utils_norec[lp] = val_lp_norec
            end

            loc_shocks = rand(dist_ζ, nl) #draw location shocks
            lp_choice = maximum(l_utils.+loc_shocks) #maximizing value
            lp = findmin(abs.(l_utils.+loc_shocks.-lp_choice))[2] #index of choice
            region_p = regions[lp]
            μ_a_R_p = μ_a_opts[region_p] #check

            lp_choice_norec = maximum(l_utils_norec.+loc_shocks) #maximizing value
            lp_norec = findmin(abs.(l_utils_norec.+loc_shocks.-lp_choice_norec))[2] #index of choice

            #update utility
            util += β * loc_shocks[lp] #location shock
            util -= (Δ_mod - Δ_2 * (h3p_val) - Δ_3 * prox_mat[l2, lp] - Δ_4 * (S-1) - Δ_5 * amenity[lp]) * (l2!=lp) * β #moving cost

            #####draw marriage/fertility outcomes
            ##marriag
            #draw
            marr_draw = rand()
            mp = 1 #non-married
            if marr_draw<prob_marr
                mp=2 #married
            end

            ##fertility
            prob_fert = prob_fert_marr
            if mp == 1 #non-married
                prob_fert = prob_fert_nonmarr
            end

            fert_draw = rand()
            fert = 0 #no kids
            if fert_draw<prob_fert
                fert = 1 #kids
            end

            #draw grandchild ability
            draw_app, app = rand(), 1 #child ability draw

            #determine child abiltiy
            tap1 = app_probs[1, region_p]
            tap2 = tap1 + app_probs[2, region_p]
            tap3 = tap2 + app_probs[3, region_p]
            tap4 = tap3 + app_probs[4, region_p]

            if draw_app>tap1 && draw_app<tap2
                app=2
            elseif draw_app>tap2 && draw_app<tap3
                app=3
            elseif draw_app>tap3 && draw_app<tap4
                app=4
            elseif draw_app>tap4
                app=5
            end

            tp = interp_t_2010_cfact(h3p, lp, app, S, R) #parent investment
            xp = interp_x_2010(h3p, lp, app, S, R) #parent goods


            #compute earnings with more flexible child marriage setup. Could technically add stuff to this, but this is already pretty grueling
            kid_marr_draw = rand()
            row_kid_marr_dist = 6*(l2-1) + 3*(S-1)
            spouse_kid_type, μ_hc_sp, σ_hc_sp = 1, 0, 0
            threshold_kid_marr_1 = kid_marr_dist[row_kid_marr_dist+1, 4]
            threshold_kid_marr_2 = threshold_kid_marr_1 + kid_marr_dist[row_kid_marr_dist+2, 4]
            if kid_marr_draw>threshold_kid_marr_1 && kid_marr_draw<threshold_kid_marr_2
                spouse_kid_type, μ_hc_sp, σ_hc_sp = 2, kid_marr_dist[row_kid_marr_dist+2, 5], kid_marr_dist[row_kid_marr_dist+2, 6]
            elseif kid_marr_draw>threshold_kid_marr_2
                spouse_kid_type, μ_hc_sp, σ_hc_sp = 3, kid_marr_dist[row_kid_marr_dist+3, 5], kid_marr_dist[row_kid_marr_dist+3, 6]
            end

            earn_sp_kid = 0
            earn_sp_kid_p4 = 0
            x_sp_kid = 0
            time_sp_kid = 0.13
            hc_draw_sp = rand(LogNormal(μ_hc_sp, σ_hc_sp))
            h3p_sp = get_index(hc_draw_sp, h_grid)
            h3p_shocks_sp = [rand(LogNormal(μ_ε3[1], σ_ε3[1])), rand(LogNormal(μ_ε3[2], σ_ε3[2]))]
            if spouse_kid_type>1 #update spousal edarnings if working
                time_sp_kid = interp_t_2010_cfact(h3p_sp, lp, app, spouse_kid_type-1, R)
                x_sp_kid = interp_x_2010(h3p_sp, lp, app, spouse_kid_type-1, R)
                earn_sp_kid = hc_draw_sp * exp(log(prim_2010.skill_prices[lp, spouse_kid_type-1]+w_mods[R])) * (1-time_sp_kid * fert)
                earn_sp_kid_p4 = exp(log(skill_prices[lp, spouse_kid_type-1]+w_mods[R])) * h3p_sp * h3p_shocks_sp[spouse_kid_type-1]
            end

            skill_price_starting = exp(log(prim_2010.skill_prices[l, S]+w_mods[R]))
            skill_price_ending = exp(log(prim_2010.skill_prices[lp, S]+w_mods[R]))
            e3p_fam = h3p_val*skill_price_ending * (1-tp*fert) + earn_sp_kid * (mp-1) #child family earnings
            e3p = h3p_val*skill_price_ending * (1-tp*fert) #compute next-period individual earnings

            ########udpate utility and consumption again -- this time a bit involved . . .#####
            #period 3 consumption, utility
            ctemp = (e3p * (1-taxes[lp]))/prices[lp] - xp
            if fert == 1
                ctemp/=1.5 #equivalence scale adjustment
            end
            con += ctemp
            util += (utility(ctemp) + (α_1 * amen_1[lp] + α_2 * amen_2[lp] + α_3 * amen_3[lp])) * (1 + θ*(mp-1)) * β #add to utility, with altruistic adjustment if married

            #realize period-4 HC
            h4p_val = h3p_val * rand(LogNormal(μ_ε3[S], σ_ε3[S]))
            h4p = get_index(h4p_val, h_grid)
            draw_college_g = rand(Normal(0, σ_η)) #college preference draw for grandkid

            e3p_fam = h3p_val*prim_2010.skill_prices[lp, S] * (1-tp*fert) + earn_sp_kid * (mp-1) #child family earnings
            #e3p_fam = h3p_val*prim_2010.skill_prices[lp, S] * (1-tp*fert) * mp #child family earnings
            e3p = h3p_val*prim_2010.skill_prices[lp, S] * (1-tp*fert) #compute next-period individual earnings
            #e3p = e3p_fam/mp #compute next-period individual earnings
            mig = (l!=lp) * 1.0

            #additional migration variables
            mig2 = (l2!=lp) * 1.0 #migrated in period 2-3

            ###update simulated data vector
            line = [fips[l] fips[lp] prim_2010.skill_prices[l, S] prim_2010.skill_prices[lp, S] m mp t e2p e3p h3p_val ap S_par S pops[l] mig type mig2 fips[lp_norec] e3_fam e3p_fam]
            push!(data_simul, line)
        end
    end

    #wrap up and export
    data_simul = vcat(data_simul...) #stick all together into one array
    data_simul_p25 = vcat(data_simul_p25...) #stick all together into one array
    pop_weights = fweights(data_simul[:,14]) #define population weights

    #compute kid and parent percentiles, with these being the final two entries of the data
    ncol = length(data_simul[1,:])

    nrows = length(data_simul[:,1]) #useful for later
    ranks = zeros(nrows, 1)
    for i = 1:nrows
        ranks[i,1] = findmin(abs.(kid_pctiles .- data_simul[i,ncol]))[2]
    end

    data_simul = hcat(data_simul, ranks)
    #data_simul_p25 = hcat(data_simul_p25, ranks[:,1])
    output_simul = data_simul
    output_simul #return simulated output
end
